/*
 * heliumatom1d.c
 *
 *  Created on: Mar 23, 2011
 *      Author: nvep5898
 */

#include "fdm1dutils.h"

main(){
	printf("inside the helium atom\n");
	gsl_matrix* H;
	/*
	 * Eigen vector/Value declaration
	 */
	gsl_vector* pot;
	gsl_vector* eig_val;
	gsl_matrix* eig_vect;
	gsl_vector* pd;
	gsl_vector* hartree_new;
	gsl_vector* hartree_old;

	initConstants();
	printConstants();

	printf("forming the potential:\n");
	pot = formPotHelium1D();


	printf("Before forming the HMatrix\n");

	/*
	 * have to loop through until the error value between
	 * the old and the new SCF/Hartree term is less.
	 *
	 */
	double error = 1;
	hartree_old = formInitialHartreePotential();
	printVector(hartree_old, "hartree old - init");

	while(error > 0.01){

		/*
		 * 1. Form HMatrix using the nuc pot and the hartree.
		 * 2. solve for eig values and vectors
		 * 3. get the probab density from the eig vectors
		 * 4. calculate new hartree from the probab density
		 * 5. calculate error func from new hartree and old hartree
		 */

		/******** STEP 1 ***********************/
//		printf("Before forming the HMatrix\n");
		H = formHeliumHMatrix(pot,hartree_old);
//		printf("Done forming the HMatrix\n");
	//	printHMatrix(H);
//		saveMatrix(H,"HmatrixHydr1d.dat");

		/******** STEP 2 ***********************/
		eig_val = gsl_vector_alloc(npoints);
		eig_vect = gsl_matrix_alloc(npoints,npoints);

		//get eigen values and eigen vectors of Hmatrix
//		printf("Solving for eigen values\n");
		solveEig(H,eig_val,eig_vect);

		/******** STEP 3 ***********************/
//		printf("Got Eig, getting probabDensity for groundlevel>> %d\n",eig_val->size);
		pd = getProbabDensityFromEigVect(eig_vect,0);

		/******** STEP 4 ***********************/
//		printf("forming the new hartree potential:\n");
		hartree_new = formHarteePotential(pd);

		/******** STEP 5 ***********************/
//		printf("calculating the error for this iteration\n");
		error = getErrorFromHartreePotentials(hartree_old,hartree_new);
		printf("this is the calculated error function:%f\n",error);
		//set the old as the new hartree
		hartree_old = hartree_new;
		getchar();
	}


	//print eigen values
	//printEig(eig_val,eig_vect);
	printf("out of the loop, making plot files\n");
	makePlotFiles(eig_val,eig_vect,"He-eigvect.dat","He-eigval.dat","He-prob.dat");

	printf("Exiting program!\n");

}

